Rigorous mean-field dynamics of lattice bosons: Quenches from 
the Mott insulator 
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Abstract. - We provide a rigorous derivation of Gutzwiller mean-field dynamics for lattice bosons, 
showing that it is exact on fully connected lattices. We apply this formalism to quenches in the 
interaction parameter from the Mott insulator to the superfluid state. Although within mean-field 
the Mott insulator is a steady state, we show that a dynamical critical interaction Ud exists, such 
that for final interaction parameter Uf > Ud the Mott insulator is exponentially unstable towards 
emerging long-range superfluid order, whereas for Uf < Ud the Mott insulating state is stable. 
We discuss the implications of this prediction for finite-dimensional systems. 
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Introduction. — Because the energy scales in ul- 
tracold quantum gases are much smaller than in solid- 
state systems, time scales are much larger. This 
provides the unique opportunity to investigate out-of- 
equilibrium many-body quantum mechanics on experi- 
mentally tractable time-scales. Moreover, these systems 
are well isolated from the environment, such that no de- 
coherence destroys the quantum correlations. Many fas- 
cinating examples have already been published, includ- 
ing collapse and revival dynamics of strongly correlated 
bosons [T], relaxation dynamics in one dimension [21 [3|, 
particle transport of fermions [4] and bosons [5] , diffusion 
dynamics of strongly correlated fermions [5], and many- 
body Landau-Zcner dynamics . 

The theoretical description of the dynamics of quan- 
tum many-body systems is very challenging, since already 
systems in thermal equilibrium require significant effort. 
However, for lattice bosons a simple, non-perturbative 
mean-field theory is available, consisting of a mean-field 
approximation in the hopping part of the Hamiltonian 
[5Hl2). Henceforth we will refer to this as Gutzwiller 
mean-field theory. This approximation can be shown to 
be exact on fully connected lattices [8] and in the limit 
of infinite dimensions [13l[T4] regardless of the strength 
of the interaction, which means that the method is fully 
non-perturbative. Gutzwiller mean-field theory offers a 
straightforward extension towards out-of-equilibrium sit- 
uations |15| , which has already been applied to many prob- 
lems: creation of a molecular condensate [15] . dynamics 



of the superfluid-insulator phase transition [TB] , transport 
in inhomogeneous systems |17| . collapse and revival os- 
cillations [T5] , atom lasers [TH] , ramp- up dynamics of the 
optical lattice p0ll21| . Bragg scattering [22], and trap dy- 
namics |23j . So far this method has been justified by semi- 
classical arguments, but no rigorous derivation is avail- 
able, making it questionable in which circumstances this 
approximation can be trusted. 

In this Letter, we show that the dynamical Gutzwiller 
equations can be rigorously derived on fully connected lat- 
tices, by using the insights recently obtained in Ref. [24] . 
This is an intriguing result, because the static Gutzwiller 
equations are also exact on the fully connected lattice. 
Moreover, static Gutzwiller mean-field theory has proven 
to be a good approximation in three spatial dimensions. 
Therefore, this results gives confidence that the dynamic 
Gutzwiller equations are also a good approximation in 
three spatial dimension. 

As an application we consider quenches from the su- 
perfluid phase towards the Mott insulator phase. Within 
Gutzwiller mean-field theory the Mott insulating state is 
trapped: the state is inert when the interaction is de- 
creased towards the superfluid regime. Here we investi- 
gate the stability of this steady state after the quench, 
by applying a small perturbation to the Mott state. We 
find that apart from the equilibrium critical interaction 
Uc = (3 + \/8)J, which separates the superfluid (for 
U < Uc) from the Mott insulator (for U > Uc), a sec- 
ond dynamical critical interaction t/^; = (3 — \/8)J exists. 
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separating the regime where the Mott state is stable (for 
final interaction Uf < Ud) from the regime where the Mott 
state is unstable with respect to emerging superfluid order 
(for Ud < Uf < Uc)- Wc can rigorously prove this when 
we use a cut-off on the particle number per site -/Vmax ^ 2 
by using the mapping to the classical system developed in 
|24| . This boundary remains fully intact when removing 
this constraint, which is shown by explicit time-evolution 
of the Gutzwiller mean-field equations. We close by dis- 
cussing the experimental implications of our findings. 

Model. — We consider the homogeneous Bose- 
Hubbard model, described by the Hamiltonian: 

U = - ^Jijb^J)] -^^rij + ^^n,(n,; - 1) (1) 

i^j i 

Here S^-^'' are the annihilation (creation) operators at site 
i and fii = blbi are the corresponding number operators. 
The particle density is controlled by the chemical poten- 
tial /i and the on-site interaction is described by the in- 
teraction parameter U. Jij contains the hopping ampli- 
tudes between the sites i and j. When applied to a finite- 
dimensional lattice, we take Jij ^ J/z when i and j are 
nearest neighbors (z denotes the number of neighbors) and 
zero otherwise. We also consider a fully connected lattice, 
on which the particles can tunnel from every site to every 
other site with equal tunnel amplitude Jij = J/V, where 
V is the number of lattice sites. 

Method. — When treated in Gutzwiller mean-field 
approximation, the Hamiltonian in Eq. ([1]) is replaced by 
a sum of single-site Hamiltonians Hmf ~ X]j ^mf with 

^MF = -J + 0*6. - \'^f) - M". + \n,(n, - 1), (2) 

together with the self-consistency relation (p = {bi). Note 
that due to the choice of the Jij the form of the mean-field 
Hamiltonian in Eq. ([2]) is independent of the lattice type. 
Since the mean-field Hamiltonian is a sum of commuting 
single-site Hamiltonians, its eigenstates are product states 
over the lattice sites: |$)mf = Yii l*^)*- The local state 
\(j>)i can be decomposed into Fock states \n)i — b] / y/ri^.\0) : 
\(f>)i = J2n=o^ '^n\n')i, where we have assumed a homoge- 
neous state, in which the coefficients c„ do not depend 
on the site index i. Wc also have explicitly included the 
cut-off in the particle number iVmax in the summation. 

The dynamical Gutzwiller equations follow from the 
Schrodinger equation with the mean-field Hamiltonian 
^ applied to the product state |$)mf: idt\^{t))MF = 
^mf|^'(0)j together with the self-consistency equation 
(j){t) ^ {^{t)\bi\^{t)) . Here and in the following we set 
ft = 1. This leads to time-dependent coefficients c„(t), 
which evolve according to the non-linear differential equa- 
tion: 



ic„{t) = -j|0(t)V^c„_i(t) +0(t)Vn+ lc„+i(t)} 
+ \^^n{n-l)- ^7ijcn{t), (3) 

together with (p{t) = X^^^^T "v/j c*_2(t)c„(t) These equa- 
tions can straightforwardly be implemented in matrix- 
form with low numerical effort, even for large values of 
^max- We exploit this in the second part of the paper, 
where we investigate quench dynamics. 

To make contact with the exact dynamics on the fully 
connected lattice, we write c„ = a„e'"", where an is the 
absolute value of c„ and q;„ its phase. We also introduce 
m„ = a^, which is the probability of having exactly n 
particles per site, or equivalently, the fraction of sites with 
exactly n particles. From Eq. ([3]) we can directly derive 
the equations of motion for Oj and rrij {Aaj = aj — aj-i): 

aj = - 1)U + ^ij 

+J ^ cos(Aq;j - Aafc) 



+ ^m.j+i{j + l)cos(Aaj+i - Aafe)|, (4) 
rrij — 2J^^ mj-mknik-ik^y^mj-ij sin(Aaj — Aak) 

k 

-^mj+i{j 4- 1) sin(AQ;j+i - Aafc)|. (5) 

On fully connected lattices the dynamics can be de- 
rived without any approximation. Here we closely fol- 
low the derivation of [53]: on fully connected lattices 
ground states must be completely symmetric when sites 
are interchanged. We can therefore build up a complete 
set of states contributing to the ground state by con- 
sidering states with fixed numbers {Mi} of sites with 
i particles, which are fully symmetrized over all sites. 
Later on we will use = Mi/V as the fraction of 
sites with i particles. Note that this definition of 
coincides with the one used before in the context of 
Gutzwiller dynamics. The states can thus be denoted 
by |Mo,...,My) and the ground state has the form 
|gs> =Y.{,ni}^Mo,Mi,---MvWn,---,Mv). Normalization 
imposes the contraint X){Afi} I^Afo,- -,Mv P = 1j whereas 
particle number conservation yields the second constraint 
Y.{Au]iY.i^'^i)\^Mo---Mv? = where N is the number 
of particles. Like in the case of the Gutzwiller mean- 
field equations, we can put a cut-off on the maximal 
number of particles per lattice site iVmax- The sim- 
plest non-trivial case is Mnax = 2. Due to the two 
constraints, the states can then be labeled by the sin- 
gle integer M2 [H]- We will make use of this simplifi- 
cation later. The symmetry constraints used to deter- 
mine the ground state also apply when studying time- 
evolution induced by a quench in the parameter U , since 
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this preserves the permutation symmetry of the Hamil- 
tonian. For general A^max, the Bose- Hubbard Hamilto- 
nian couples states | . . . , Afj_i, Mj, . . . , Mk-i,Mk, ■ ■ ■) to 
, M,_i -h 1, M, - 1, . . . , Mt-_i - 1, Mfc -I- 1, . . .) with 

r jj are the diagonal elc- 

. , mj_i, TTij, TTij + i, . . .) to 

.). The Schrodingcr 
expressed in terms 



amplitude Wjk 
ments and VFjj+i couples 



Note that W,. 



, mj_i-|-l,TOj -2, . . . , nij+i + l, . 



■Mv 



equation for the amplitudes Amq ■ 
of the Wjk therefore yields: 



iAMo---M,^iMyMk-iMk---Mv - (6) 
W"jfc^Mo---(A/,_i-H)(M,-l)---(j\4_i-l)(j\4-H)---Afv 

jk 

We now consider the thermodynamic limit V ^ oo. In 
this limit the transition amplitudes Wki assume a simple 
form: 



W. 



jk 



U 



J^MjMj-iMkMk 




Ojk 



(7) 



such that Wjk ~ Wkj ■ Tm'ning now to the {wj} variables, 
which become continuous in this limit, the Schrodingcr 
equation turns into: 



iA{{mj}) = Wjk cos{pj 
Jk 



Pj^i -Pk-i +Pk)A{{mj}) 



where we have introduced for each rrij an associated mo- 
mentum pj = —idmjy in exactly the same manner as 
was derived in |24j for A^max = 2. In the thermodynamic 
limit nij and Pj commute and we obtain classical equa- 
tions of motion nij = dH/dpj and Pj = —dH/dnij for the 
classical Hamiltonian 



H 



U 



"ijJ - (8) 

J j 
J y^ y/mjmj^imknik-i\/jk cos{Apj - Apk) 

It is now straightforward to see that these classical equa- 
tions of motion exactly coincide with the Gutzwiller equa- 
tions of motion in Eqs. (jH [S]) when we identify the phase 
aj with the momentum pj . 

This constitutes the proof that Gutzwiller dynamics is 
exact on fully connected lattices. Therefore it is a con- 
trolled mean-field method with the same range of validity 
as the static Gutzwiller equations which are also exact on 
fully connected lattices [5]. 

The fact that the static and dynamic Gutzwiller equa- 
tions are both exact on fully connected lattices, i.e. in 
infinite dimensions, gives support for the idea that the 
accuracy of this mean-field method in finite dimensional 
lattice is comparable for both cases. The accuracy of the 
static mean-field theory has been extensively tested, be- 
cause the results can be compared with numerically exact 
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Fig. 1: Classical paths for A'^max = 2 with H = for various 
values of Uf. 



results from Quantum Monte Carlo simulations PSII^ . 
One can for instance compare the prediction for the crit- 
ical interaction {U/J)c for the supcrfluid-Mott insulat- 
ing transition. This yields reasonably good quantitative 
agreement in three dimensions. In two dimensions the 
agreement is merely qualitative. Also for inhomogeneous 
systems Gutzwiller mean-field predictions have been com- 
pared with Quantum Monte Carlo data, yielding a similar 
agreement as for the homogeneous case [57] ■ This pro- 
gressively good accuracy of the mean-field theory when 
the dimension is increased, can be understood by viewing 
the static Gutzwiller mean-field equations as the leading 
terms in a \/z expansion [28l[29]. Including the second 
order terms yields Bosonic Dynamical Mean-Field The- 
ory p8H32] . for which almost exact agreement with the 
Quantum Monte Carlo results has been obtained [3T1I32] . 
Several efforts have already been made to go beyond the 
mean-field theory for the out-of-cquilibrium situation as 
well ^34^ 

It is worth noting that although the mean-field theory 
turns out to be a good approximation in sufficiently high 
dimensions, the ground state description does not include 
any entanglement between different lattice site: in mean- 
field the states are perfect product states. 

The derivation presented here is clearly independent of 
the strength of the interaction parameter and hence valid 
both in weak and strong coupling. The equations are 
moreover valid for all types of dynamics: since no adi- 
abatic assumptions have been made, the Gutzwiller dy- 
namics can be applied for slow, quasi-adiabatic dynamics, 
as well as for fast dynamics. An example of the latter is 
the quench dynamics we turn to now. 

Quenches from the Mott insulator. We study 
the situation in which after preparing the ground state for 
some value oiU = Ui, the interaction is suddenly changed 
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toU — Uf at time t = and the subsequent time-evolution 
of the double occupancy 7712 (t) and the condensate fraction 

N^{t) = Y.t^M^i^/'^'^ = l'?^'^^)!^ studied. We study the 
case that the particle density n = N/L = 1. As is well 
known, the equilibrium properties are then that within 
Gutzwiller the system is superfluid for U < Uc = {'i + 
V8)J and Mott insulating for U > Uc- This last state is 
particularly featureless withing Gutzwiller: it is a perfect 
product state with exactly one particle per site. 

In previous work }24) quenches within the superfluid 
phase were studied. A dynamical transition was found 
between small and large quantum quenches. 

Here we study quenches from the Mott insulating state, 
i.e. we choose Ut > Uc [34H36] . This might seem prob- 
lematic, because the Mott state is completely inert within 
mean-field theory: when the system is in the Mott insu- 
lating state, the system will remain in this state for all 
Uf. We overcome this problem here, by analyzing the 
stability of the Mott state after a quench: a small per- 
turbation is added to the perfect Mott state, which can 
drive the system away from the insulator and thus prove 
the insulator unstable. We note that within mean-field 
this perturbation might seem artificial; however, quantum 
fluctuations beyond mean-field, which are always present 
in finite-dimensional systems, naturally provide this per- 
turbation. 

We now first analyze the situation for A'^max = 2. In 
this case, the equations of motion become very simple, 
because of the conservation of particle number and norm. 
Using this, one obtains coupled classical equations of mo- 
tion for (m2,P2) with H = Um2 — J{1 — 2m2)m2(3 4- 
V8cos(]32)) El. Moreover A^cW = m2(i)(l - 2m2(i))(3 + 
V8cosp2{t)), and hence it is only nonzero if 1712 is nonzero 
and proportional to m2{t) for small m2{t). Therefore we 
focus on the time evolution of m,2{t) in the following. 

For the Mott insulating state 7712 = and hence the 
energy H = 0. Since energy is conserved, the system has 
to evolve along paths with H = 0. These are shown in Fig. 
[T]for various values of Uf. The path with m2 = is always 
a solution, and when initially 7712(0) = 0, this will remain 
true for all later times. (This is true in the thermodynamic 
limit; for finite systems on the fully connected lattice the 
system can escape from this state on a timescale diverging 
with the system size as \og{V) [M].*) 

We now want to investigate the stability of this solution: 
what happens if 777-2(0) is slightly positive? To answer 
this question we look at the other branch with H = 
which is characterized by 7772 > 0. Now a fundamental 
difference appears between Uf < Ud and Uf > Ud (where 
Ud = {3 — \^)J): in the former case this branch is fully 
disconnected from the 7712 = curve. This implies that 
paths with initial conditions with 7772 <^ 1 remain infinites- 
imally close to 777,2 = and the Mott state is dynamically 
stable under small perturbations. An important direct 
consequence of this is that the system does not thermal- 
ize: because the system is completely static, there is no 
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Fig. 2: Time evolution of the double occupancy m2{t) starting 
from the Mott insulator witli a tiny perturbation quenched 
towards Uf > Ud {Uf = J) and Uf < Ud {Uf = 0.15J). Shown 
are results from Gutzwiller dynamics (G.) with different A^max, 
and classical dynamics (CI.) for Nrmai = 2. 

evolution towards a state that could be characterized by 
a fictitious temperature. 

In contrast, when Uf > Ud, the two branches cross, 
such that in this case solutions with 7772(0) > will fol- 
low a curve on which 7712 (t) grows macroscopically large. 
Hence the Mott state is unstable for this case. We verify 
this numerically by time-evolving the classical equations of 
motion and the Gutzwiller equations for initial states that 
arc slightly perturbed away from the pure Mott insulator. 
This fully confirms this picture: for Uf > Ud the double 
occupancy grows exponentially, whereas for Uf < Ud it re- 
mains small and oscillates. Note that these oscillations are 
fully within the regime of an infinitesimal 777-2 and hence a 
stable Mott insulator. They correspond to the time evolu- 
tion along the path infinitesimally close to the line 7772 = 0, 
which still involves dynamics oi p2, refiecting itself in os- 
cillations in 7712 as well. 

This is shown in Fig. [21 in which also results with higher 
Amax are shown. We observe that the initial exponential 
growth is independent of A^max^ only after the superfluid 
order parameter has grown macroscopically large, differ- 
ences due to different values of A^max are visible. We can 
fit the double occupancy with an exponential fit over many 
orders of magnitude: 7712 (t) cx; e'*'*''. The exponent A from 
this fit is shown in Fig. |21 It can be derived analyti- 
cally by investigating the classical equations of motion for 
A^max = 2. The time-evolution starts with 7772 = P2 = 0. 
First the system evolves along the 7772 = 0-line, until 
P2 = at cosp ^ {Uf — 3J)VSiJ. Then the curve sharply 
bends and moves with constant p2 driven by the linearized 
equation of motion 

7772 = V8Jm2 sinp2 = 7712^^ {Uf — Ud){Uc — Uf). 

This yields the observed exponential increase, with expo- 
nent A = \/{Uf — Ud){Uc — Uf), which exactly fits the nu- 



p-4 



Rigorous mean-field dynamics of lattice bosons: Quenches from the Mott insulator 



2.5 
2 

1.5 
1 

0.5 




Numerical 
Analytical 



U,/J 



Fig. 3: Exponent A of the exponential growth of m2{t) oc 
e^'"' as a function of Uf, both obtained by numerical fitting 
(squares) and the analytic expression (line). 



merical fits in Fig. |31 We observe a remarkable symmetry 
around Uf = 3 J in this picture. The exponent vanishes 
like a square-root when Uf approaches Uc and Ud, like one 
expects from mean-field dynamics. It is worth noting that 
the critical point does not shift in this case, as generally 
expected when studying quench dynamics in dimensions 
d> 2 [37] . This is because the mean-field Mott insulating 
state we start from is a perfect product state, in which the 
correlation length exactly vanishes. 

Conclusions and experimental implications. 

We have shown that Gutzwiller mean-field theory is 
exact on fully connected lattices, i.e. in infinite dimen- 
sions, irrespective of the interaction strength or the ques- 
tion whether the dynamics is slow or fast. We applied this 
mean-field formalism to quenches in the ?7/J-ratio from 
the Mott insulator to the superfiuid side of the phase dia- 
gram. Although the pure Mott insulator is always a steady 
state, we found that infinitesimal fiuctuations make this 
state exponentially unstable towards emerging long-range 
superfiuid order for Uf > Ud- 

Whereas this is a rigorous statement on the fully con- 
nected lattice, we expect this dynamical transition to be 
present on lattices with finite connectivity as well. Espe- 
cially for three-dimensional lattices, where quantum fiuc- 
tuations only lead to small quantitative corrections to the 
static Gutzwiller predictions, we expect the Gutzwiller dy- 
namics to be a good approximation and the dynamical 
transition to be present. A second argument for this comes 
from the consideration of the special point Uf — 0. In this 
case the Hamiltonian driving the time-evolution consists 
only of the hopping part, which commutes with operators 
of the form b^^bi measuring long range phase coherence. 
Hence it is a rigorous statement that for J7/ = no su- 
perfiuid coherence can build up and the Mott insulating 
state is stable for this point. Whereas Ud will be shifted 
away from the mean-field value in finite dimensional lat- 
tices (like also the value of Uc is changed) , this argument 



thus implies that Ud > 0- However, even the case Ud = 
would be interesting, because in principle it is also possible 
to quench to attractive interactions. 

Another important aspect that changes in finite dimen- 
sions, is that the Mott insulating state is no longer feature- 
less: virtual particle hole-pairs lower the energy propor- 
tional to J^/C/i. This naturally provides the perturbation 
from the perfect product state that is needed to trigger 
the exponential build-up of superfiuid order, but it also 
means that the exact phase boundary can only be probed 
for initial Ui/J^ 1. 

The existence of a dynamical transition can experimen- 
tally be investigated by using ultracold atoms in optical 
lattices, in which the ratio U /J can easily be quenched by 
modifying the depth of the optical lattice or by means of 
Feslibach resonances. In this latter case it is also possi- 
ble to quench to attractive interactions. The occurrence 
of long-range superfiuid order is directly visible in time-of- 
fiight experiments, in which the phase coherence is probed, 
and also the amount of double occupancies can be directly 
observed |38j . 

However, the presence of a trapping potential, leading to 
an inhomogeneous system, complicates the situation. This 
is because of two reasons: after a quench in the U/ J-ratio 
also particle transport sets in, such that the local ener- 
gies and local particle numbers are no longer conserved. 
Secondly, if the Mott insulator is surrounded by a super- 
fiuid, phase coherence is build up from the outside. The 
sharp dynamical transition is therefore best visible in a 
fiat potential, as realized in [6]. 

This work was supported by the Nederlandse Organ- 
isatie voor Wetenschappelijk Onderzoek (NWO). 
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